Non-Gaussian energy landscape of a simple model for strong network-forming liquids: 
Accurate evaluation of the configurational entropy. 
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We present a numerical study of the statistical properties of the potential energy landscape of 
a simple model for strong network-forming liquids. The model is a system of spherical particles 
interacting through a square well potential, with an additional constraint that limits the maximum 
number of bonds, iV max , per particle. Extensive simulations have been carried out as a function of 
temperature, packing fraction, and iV max . The dynamics of this model are characterized by Arrhe- 
nius temperature dependence of the transport coefficients and by nearly exponential relaxation of 
dynamic correlators, i.e. features defining strong glass-forming liquids. This model has two impor- 
tant features: (i) landscape basins can be associated with bonding patterns; (ii) the configurational 
volume of the basin can be evaluated in a formally exact way, and numerically with arbitrary pre- 
cision. These features allow us to evaluate the number of different topologies the bonding pattern 
can adopt. We find that the number of fully bonded configurations, i.e. configurations in which all 
particles are bonded to iV max neighbors, is extensive, suggesting that the configurational entropy of 
the low temperature fluid is finite. We also evaluate the energy dependence of the configurational 
entropy close to the fully bonded state, and show that it follows a logarithmic functional form, differ- 
ently from the quadratic dependence characterizing fragile liquids. We suggest that the presence of 
a discrete energy scale, provided by the particle bonds, and the intrinsic degeneracy of fully bonded 
disordered networks differentiates strong from fragile behavior. 
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I. INTRODUCTION 

When a liquid is fastly supercooled into a metastable 
state under the melting point, its structural relaxation 
time, t, increases over 13 orders of magnitude with de- 
creasing temperature, T. Below some given temperature, 
equilibration is not possible within laboratory time scales 
and the system becomes a glass Q, . The glass tran- 
sition temperature, T g , is operationally defined as that 
where r = 100 seconds, or the viscosity rj = 10 13 poise. 

Angell has introduced a useful classification scheme Q 
for glass-forming liquids. According to the definition of 
kinetic fragility a liquid is classified as 'strong' or 'frag- 
ile' depending on how fast its relaxation time increases 
when approaching T g . Liquids that show a weak depen- 
dence, well described by an Arrhenius law r oc exp(^4/T), 
with A a temperature independent quantity are classi- 
fied as strong. Strong liquids form open network struc- 
tures that do not undergo strong structural changes with 
decreasing temperature. On the contrary the dynam- 
ics of fragile liquids, as many polymeric or low molec- 
ular weight organic liquids, where interactions show a 
less directional character, are more sensitive to temper- 
ature changes, and relaxation times show strong devia- 
tions from Arrhenius behavior 0, Q . Several empirical 
functions have been proposed for the T dependence of 
r in fragile liquids, the Vogel-Tammann-Fulcher (VTF) 



law, t oc exp[A/ (T — Tq)], having gained more acceptance 
|(|. In this equation To is the VTF temperature. 

Kauzmann noted |?J that, when extrapolating to low 
temperature the experimental T dependence of the con- 
figurational entropy, Sconf, the latter became zero at 
a certain temperature Xk ('Kauzmann temperature') 
somewhere below T g . Experiments Q, § often provide 
the result Tk ~ To. Given the Arrhenius character of 
strong liquids, this comparison would also suggest that 
Tk is zero for these systems. 

However, it must be stressed that the values of Tk 
and To are the result of an extrapolation of experimen- 
tal data, which in principle is not necessarily correct. In 
particular, an extrapolation below Tk would lead to an 
'entropy catastrophe': a disordered liquid state with less 
entropy than the ordered crystal. In practice, the equi- 
librium liquid state at the putative Tk is never reached 
in experiments, because the liquid falls out of equilibrium 
at T g > Tk- The fate of the configurational entropy in an 
ideal situation where arbitrarely long equilibration time 
scales could be accessed is one of the key (and contro- 
versial) questions associated to the glass transition prob- 
lem. One solution states that crystallization is unavoid- 
able when approaching Tk in equilibrium 0, flo| . It has 
also been proposed that Sconf changes its functional form 
below T g , remaining always positive |ll| . Another solu- 
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tion is that ffronf reaches zero at 7k and remains constant 
below it ElGlElEi 

Some insight into this latter question, in the physical 
origin of the fragility, and in general in the relation be- 
tween dynamic and thermodynamic properties of glass- 
forming liquids, can be obtained by investi gating th e po- 
tential energy landscape (PEL) EH El El ElM E3, 
i.e., the topology of the potential energy of the liquid 
U = U(r ). According to the inherent structure (IS) for- 
malism introduced by Stillinger and Weber E3 > the PEL 
is partitioned into basins of attraction around the local 
minima of U. These minima are commonly known as the 
'inherent structures'. The free energy is obtained as a 
sum of a 'configurational' contribution, resulting from 
the distribution and multiplicity of the different IS's, 
and another 'vibrational' contribution, resulting from the 
configurational volume available within the basin around 
each individual IS. The introduction of the IS formula- 
tion has motivated a great theoretical and computational 
effort in order to understand the connection between the 
statistical properties of the PEL and th e dy namic behav- 
ior of supercooled liquids jl|2i|^|2l[ H |2 i | ^ |2 l |3 ^ 

Enl 0130013 III 

and nowadays has become a key methodology in the field 
of the glass transition. 

From a series of numerical investigations in models of 
fragile liquids, it is well established that for such systems 
the distribution of inherent structures is well described 
by a Gaussian function, at least in the energy range that 
can be probed within the equilibration times permitted 
by computational resources [2^, 0, 0] . It can be for- 
mally proved E! 0] that for a Gaussian landscape the 
energy of the average visited IS depends linearly on T , 
a result that has been verified in several numerical stud- 
ies of fragile liquids 0, 0, 0, 0, 0] . Recent studies of 
the atomistic BKS model for silica |3ll lijj , the archetype 
of strong liquid behavior, have shown instead that devia- 
tions from Gaussianity take place in the low energy range 
of the landscape, and indeed, at low temperature, the av- 
erage inherent structure progressively deviates from lin- 
earity in T^ 1 . The existence of a lower energy cut-off and 
a discrete energy scale, as it would be expected for a con- 
nected network of bonds, has been proposed as the origin 
of such deviations from Gaussianity |43j ] . These investi- 
gations also suggest that S^onf cannot be extrapolated to 
zero at a finite T, and as a consequence, strong liquids 
would not show a finite 7k. However, long equilibra- 
tion times prevent the determination of the lowest-energy 
state and its degeneracy, and therefore an unambiguous 
confirmation of this result. 

Recently, we have proposed a minimal model which 
we believe capture the essence of the prototype strong 
liquid behavior. In this model particles interact 

via a spherical square-well potential with an additional 
constraint on the maximum number of bonded neigh- 
bors, A max , per particle. The lowest-energy state is the 



fully bonded network and its energy is thus unambigu- 
ously known. Within the equilibration times permitted 
by present day numerical resources, configurations with 
more than 98 per cent of the bonds formed can be prop- 
erly thermalized. As a result, no extrapolations are re- 
quired to determine the low T behavior. This model is 
particularly indicated for studying the statistical prop- 
erties of the landscape. In analogy with the Stillinger- 
Weber formalism, we propose to partition the configura- 
tion space in basins which, for the present case, can be 
associated with bonding patterns. A precise definition 
of the volume in configuration space associated to each 
bonding pattern can be provided, since a basin is charac- 
terized by a flat surface with energy proportional to the 
number of bonds. Crossing between different basins can 
be associated to bond-breaking or bond- forming events. 
In contrast to other systems previously investigated, the 
vibrational contribution of the PEL can be expressed in 
a formally exact way. No approximation for the shape 
of the basins is requested. As we discuss in the follow- 
ing, the precision of the evaluation of the basin volume, 
only limited by the numerical accuracy of statistical av- 
erages, allows us to evaluate, with the same precision, 
the configurational entropy. 

In this article we report the study of the statistical 
properties of the A max model for the cases A max = 3, 4 
and 5, and for a large range of packing fractions <f>, ex- 
tending the results limited to a fixed <p and N m xx = 4 
previously reported in a short communication 48] . The 
range of 4> values here investigated are (0.2-0.35) for 
A max = 3, (0.3-0.5) for N max = 4, and (j> = 0.35 for 
A max = 5. The article is organized as follows. In Sec- 
tion II we introduce the model and provide computa- 
tional details. In Section III we briefly summarize the IS 
formalism and apply it to the present model. Dynamic 
and energy landscape features are shown and discussed 
in Section IV. Conclusions are given in Section V. 

II. MODEL AND EVENT-DRIVEN 
MOLECULAR DYNAMICS 

The model we investigate is similar in spirit to one pre- 
viously introduced by Speedy and Debenedetti 0]. In 
the present model, particles interact through a spherical 
square-well potential with a constraint on the maximum 
number of bonds each particle can form with neighboring 
ones. Namely, the interaction between any two particles 
i and j that each have less than A max bonds to other 
particles, is given by a spherical square- well potential of 
width A and depth — u Q : 

{co r < <7 

-Uq a <r <a + A (1) 
r >a + A, 

with r the distance between i and j. When a < r < 
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a + A, particles i and j form a bond, unless at least 
one of them is already bonded to other A max particles. 
If this is the case V^{r) is simply a hard sphere (HS) 
interaction: 

W0 = {~ T (2) 
10 r > a. 

In the original model introduced by Speedy and 
Debenedetti |49j an angular constraint was imposed by 
avoiding three particles bonding loops. In an effort to 
grasp the basic structural ingredients originating strong 
liquid behavior, we have not implemented this additional 
constraint. The potential given by Eqs. |Q|2), despite its 
apparent simplicity, is not pairwise additive, since at any 
instant, the interaction between two given particles does 
not only depends on r, but also on the number of particles 
bonded to them (if a < r < a + A). Hence, to propagate 
the system it is requested not only the coordinates and 
velocities, but also the list of bonded interactions. 

We also note that the model is not deterministic. Con- 
sider a configuration in which particle i is surrounded 
by more than N max other particles within a distance 
a < r < a + A. Of course, only A max of these neighbors 
are bonded to i, i.e. feel the square-well interaction. If 
one of the bonded neighbors moves out of the square- well 
interaction range (i.e. to a distance r > a + A), a bond- 
breaking process occurs. According to Eq. QJ, a new 
bond can be formed with one of the other non-bonded 
particles whose position is in the range a < r < a + A. 
If several candidates are available — where a candidate is 
defined as a particle whose distance from i is in the range 
a < r < g + A and which is engaged in less than N max 
bonds to other distinct particles — then one of them is 
randomly selected to form the new bond with the parti- 
cle i. Of course, in each bond-breaking/formation process 
the velocities of the two involved particles are changed to 
conserve energy and momentum (see also below). 

Despite these complications, the model given by Eqs. 
Q |2l can be considered among the simplest ones for 
simulating clustering and network formation in fluids 
[50L l5l| . It does not require three-body angular forces 
or non-spherical interactions. The penalty for retaining 
spherical symmetry is the complete absence of geometri- 
cal constraints between the bonds. 

The maximum number of bonds per particle is con- 
trolled by tunning N max . For a system of N particles, the 
lowest-energy state, corresponding to the fully bonded 
network, has an energy = —NN max uo/2. If Abb is 
the number of broken bonds for a given configuration of 
the network, the energy of that configuration is given by 
E = E [h + N hh u . 

The model parameters have been set to A/(<7 + A) = 
0.03, Uo = 1 and a = 1. In the following, entropy, S, will 
be measured in units of fcs- Setting fee = 1, potential 
energy, E, and temperature, T, are measured in units of 



uq. Distances are measured in units of a. Diffusion con- 
stants and viscosities are respectively measured in units 
of a^a/m) 1 / 2 and (mu ) 1 / 2 <r~ 2 . We have simulated a 
system of TV = 10000 particles of equal mass m = 1, im- 
plementing periodic boundary conditions in a cubic cell 
of length Lbox- We have evaluated the PEL properties for 
several values of A max , in a wide range of T, and pack- 
ing fraction <j) = ttN/6L^ ox . The system does not exhibit 
phase separation for the state points here investigated 

HIE! El. 

Dynamic properties are calculated by molecular dy- 
namics simulations. We use a standard event-driven algo- 
rithm E3 for particles interacting via discontinuous step 
potentials. The algorithm calculates, from the particle 
positions and velocities at a given instant to , the instants 
t C oll and positions for all possible collisions between dis- 
tinct pairs, and selects the one which occurs at the small- 
est t co n- Then the system is propagated for a time t co \\—t 
until the collision occurs. Between collisions, particles 
move along straight lines with constant velocities. Colli- 
sions can take place at relative distance r = a (in which 
case velocities of colliding particles are reversed according 
to hard-sphere rules) and, only for bonded interactions, 
at distance a + A (in which case velocities are changed in 
such a way to conserve energy and momentum). Start- 
ing configurations are selected from previously generated 
hard-sphere configurations with hard-sphere radius cr+A. 
In this way, we start always from a configuration where 
no bonds are present. After thermalization at high T 
using the potential defined in Eqs. Q |2J), the system 
is quenched and equilibrated at the requested tempera- 
ture. Equilibration is achieved when energy and pressure 
show no drift, and particles have diffused, in average, sev- 
eral diameters. We also confirm that dynamic correlators 
and mean squared displacements show no aging, i.e., no 
time shift when being evaluated starting from different 
time origins. Once the system is equilibrated, a con- 
stant energy run is performed for production of configu- 
rations, from which diffusivities and dynamic correlators 
are computed. Statistical averages are performed over 
typically 50-100 independent samples. Standard Monte 
Carlo (MC) simulations are carried out for the calcula- 
tion of the vibrational contribution of the PEL (see Sec- 
tion III), using previously equilibrated configurations. 

III. INHERENT STRUCTURE FORMALISM 

If U = U(r N ) is the potential energy of a system 
of N identical particles of mass to, the partition func- 
tion is given by the product Z = Z 1S Z CX , where Z lg 
is the purely kinetic ideal gas contribution, and Z ox is 
the excess contribution resulting from the interaction po- 
tential. The ideal gas contribution to the free energy, 
F is = -ksTlnZte, is given by: 

/3F is /N = ln( N/V ) + 3 In A - 1 , (3) 
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where V is the system volume, (3 = (fee? 1 )" 1 , and A = 
hi^Trmk^T)^ 1 ^ 2 is the de Broglie wavelength, with h the 
Planck constant. The excess contribution to the partition 
function is: 



Z cx = j dr N exp[-pU(r N )]. 



(4) 



The IS formalism introduced by Stillinger and Weber 
(i"t| . provides a method for calculating Z cx . According 
to Stillinger and Weber the configuration space is par- 
titioned into basins of attraction of the local minima of 
the PEL, the so-called inherent structures. If £l(Eis) is 
the degeneracy of a given IS, i.e., the number of local 
minima of U(r N ) with potential energy Eis, the configu- 
rational entropy is defined as S^onf (^is) = ln[fi(.Eis)]. 
To properly evaluate Z ex , besides the information of the 
number of distinct IS's, it is necessary to evaluate the 
partition function Z°f h (T, Eis) constrained in the volume 
of each of these basins. The quantity Z°f h (T, Eis) pro- 
vides a measure of the configurational volume explored 
at temperature T by the system constrained in the basin 
of depth £"is. Averaging over all basins with the same 
depth Eis one obtains |17| : 



Z vib( T i E Is) 



1 



E 



basms(£is) 



dY " e -PiU-E ls) 



(5) 



where the notation "basins(i?is)" recalls that the sum 
is performed over all the basins whose potential energy 
minimum is Eis, and each integral runs over the config- 
urational volume associated to the corresponding basin. 
It is convenient to express Z°f h (T, Eis) as: 



(6) 



to stress that the 'free energy', f°* h (T, Eis ) , commonly 
referred as the vibrational PEL contribution, accounts 
for the vibrational properties of the system constrained 
in a typical basin of depth Eis- 

Within the above framework, the excess partition func- 
tion is given by: 



Z** = ^e 
is 



-(3[E ls -TS conl (E IS )+f%b( T > E is)] 



(7) 



In the thermodynamic limit the excess free energy, F c * — 
— fceTlnZ", can be evaluated as: 

F c X(T) = E[T) _ TSconi{E{T)) + f e? h (T,E(T)), (8) 

where E(T) is the average IS potential energy at temper- 
ature T. Therefore, the excess free energy is the sum of 
three contributions. The fist term in Eq. ||SJ accounts for 
the average value of the visited local minima of the poten- 
tial energy. The second term is related to the degeneracy 
of the typical minimum. The third term accounts for the 
configurational volume of the typical visited basin. 



The IS formalism has been applied in the past to sev- 
eral numerical studies of models of liquids. Indeed, sim- 
ulations offer a convenient way to evaluate F CX (T), E(T) 
and f°f h (T, E(T)), and to derive, by appropriate sub- 
structions, the configurational entropy. The only approx- 
imation performed in these studies refers to the vibra- 
tional free energy, which is usually calculated under the 
harmonic approximation — by solving the eigenfrequen- 
cies of the Hessian matrix evaluated at the IS — or, in 
the best cases, including anharmonic contributions under 
the strong assumption 56( that these corrections do not 
depend on the value of Eis- 

In the case of the N max model, the evaluation of the 
free energy in the Stillinger- Weber formalism is straight- 
forward. The partition function can be formally written 
as a sum over all distinct bonding patterns — i.e. over all 
configurations which cannot be transformed by deforma- 
tion to each other without breaking/forming bonds — in 
a way that is formally analogous to the IS approach once 
one identifies a bonding pattern with an IS basin. Differ- 
ently from the standard IS approach, the present specific 
step-wise potential does not require a minimization pro- 
cedure to locate the local minimum. Each bonding pat- 
tern can be associated to a basin characterized by a flat 
surface with energy proportional to the number of bonds. 
Crossing between different basins requires bond-breaking 
or bond- forming events. While in the IS formalism the 
partition function is associated to local minima, in the 
N max model Z cx is evaluated expanding around all dis- 
tinct bonding patterns. The flat surface of the basin and 
the clear-cut basin boundaries make it possible to evalu- 
ate the vibrational contribution of the PEL in a formally 
exact way, in contrast to other systems previously inves- 
tigated. No approximation for the shape of the basins 
is requested. In the A^ max model, the excess vibrational 
free energy is purely entropic. 

To calculate we make use of the Perturbed Hamil- 
tonian approach introduced by Frenkel and Ladd |57ll58| . 
which provides an exact analytical formulation for the 
excess free energy of a given system by integration from 
a reference Einstein crystal. In brief, to calculate the 
free energy of a system defined by a Hamiltonian H, 
one can add a harmonic perturbation, H\(r N ; A max ) = 

Amax J2iLi( r i — r i) 2 > around a disordered configuration 
r wo = ( r o _ it can b e demonstrated that the ex- 

cess free energies of the perturbed, F cx (T;A max ), and 
unperturbed, F ex (T; A = 0), systems are related as 

nim Hi: 

^ cx (T;A = 0)=F cx (T;A max ) 

ln[A max ] N 

A<5>-r°) 2 Mln[A]. (9) 

-°° f=i 

Brackets denote ensemble average for fixed A. Due to the 
presence of the harmonic perturbation, particles in the 
perturbed system (H + H\) remain constrained around 
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As a consequence, {J2iLi( r i ~~ r ?) 2 )A is finite. For the Carnahan-Starling formula [60j: 



a sufficiently large value of A max (so that H is negligi- 
ble as compared to H\), the perturbed system behaves 
like an Einstein crystal, i.e., like a system of 3N inde- 
pendent harmonic oscillators with elastic constant A max . 
The excess free energy for the latter system is given by: 



/?F ex (A max )/iV = (3E/N 



3 / nk B T 

2 A 

+l-bx(N/V). (10) 

In this expression E is the energy of the system in r N0 . 
If the condition \ max (E?=i(*i -r?) 2 ) Amal = 3Nk B T/2 is 
fulfilled, the harmonic limit is recovered and will be also 
valid for A| nax > A max . Hence, A max can be taken as an 
upper cut-off for the integration in A, since selecting a 
larger value A^ ax leads to trivial cancelation in Eqs. 

EDJ. 

Next we discuss how the Perturbed Hamiltonian ap- 
proach can be used for evaluating the excess vibrational 
free energy of a typical bonding pattern in the 
N max model. We start from an arbitrary equilibrium 
configuration r^ at T, and with its associated bond- 
ing pattern. We then calculate, via MC, the quantity 
(J2iLi( r i ~ r i > ) 2 )A, imposing the constraint that bonds 
can neither break nor reform. In the limit A — > the 
system samples the volume in configuration space asso- 
ciated to the selected bonding pattern. Hence, we can 
identify f°? h (T,E(T)) with F ox (T;A = 0). By retain- 
ing the constraint on the bonding pattern, we simulate 
the perturbed system for several values of A to properly 
evaluate the integrand of Eq. J5J and estimate the value 
of f°fo(T,E(T)). To improve statistics, an average over 
several starting configurations r* is performed. 

Within the above framework, the corresponding ex- 
pression for the excess vibrational entropy, S°f h , is: 



OCX o 



Ttk B T 

Xh>, 2 ~ V A max 

Mn[A max ] N 



" 1 + 111 1 7 



A<^( ri -r°) 2 ) A dln[A]. 



(11) 



We stress that Eq. I|ll|l is an exact relation for S*f h 
for the present N max model. Hence, the precision in the 
evaluation of this quantity does not depend on any ap- 
proximation, but only on the numerical accuracy of the 
MC calculation and the statistical average. 

The total excess entropy S°* t (T) can also be calculated 
with arbitrary precision by thermodynamic integration 
from a reference state at T = T Ie f where S°* t is already 
known. One possible choice is to select as reference point 
the ideal gas and integrate along a path which does not 
cross any phase boundary or, in the present case, to inte- 
grate from very high temperature. Indeed, at sufficiently 
high T rc f, the model is equivalent to a hard-sphere sys- 
tem, whose free energy is well known. An accurate es- 
timate of the hard-sphere excess entropy is provided by 



^HS 



0(30 - 4) 
(1 - 0) 2 ' 



(12) 



The excess entropy at finite T can then be obtained by 
integrating along a constant volume V path as 



ccx 
D HS 



T rof 



dE 
df 



dT 



(13) 



From the two accurate evaluations of S^t 



[Eq. PJl ] 

and S** b [Eq. ifTTjl] an accurate estimate of the con- 
figurational entropy iS C onf can be obtained as S con [ = 
&tot ~ ^vib- The resulting S^onf can be related to T or, 
parametrically, to E(T) (i.e., to the number of bonds). In 
the rest of the article, all the thermodynamic functions \P 
will be expressed as "quantity per particle" , but for sim- 
plicity of notation we will drop the factor 1/N. Hence, 
in the following "iff" will be understood as "4" /N" . 

IV. RESULTS AND DISCUSSION 

a) Dynamics 

In this Section we show that the dynamics in the JV max 
model meet the criteria defining strong liquids. The key 
feature is the Arrhcnius behavior of the transport coeffi- 
cients. Fig. H shows the T dependence of the diffusivity 
D and the viscosity r\ for different values of A^ max and 
(j). The diffusivity is calculated as the long time limit of 
(J2iLi[ r i(.t) ~ r i(0)] 2 )/6A^. The viscosity rj is determined 
as: 



'i=wb^> ((, -' 4(0)12> ' 



(14) 



where A{t) = m £^ . 



as explained in [6l|. Greek 



symbols denote the x, y, z coordinates of the position, r i: 
of particle %. An average is done over all the permuta- 
tions with a 7^ (3. As shown in Fig. ^ at low temper- 
ature both quantities display Arrhenius behavior. The 
activation energies for D and rj are approximately uo, 
suggesting that all bonds break and reform essentially in 
an independent way. This is consistent with the absence 
of angular constrains in this model. Despite the simu- 
lations are limited in time by computational resources, 
the fact that the Arrhenius form covers more than three 
orders of magnitude and the observed value of the acti- 
vation energy strongly suggest that this functional form 
will be retained at lower T. We also find (see Fig. 
that the Stokes-Einstein relation j^l, Drj/T = (37rcr) -1 , 
is fulfilled essentially at all temperatures, independently 
from the N mELX value. 

Long time decays of dynamic correlators in su- 
percooled states are usually well described by the 
phenomenological Kohlrausch- Williams- Watts (KWW) 
function exp[— (t/r)' 3 ], where r is the corresponding re- 
laxation time and p is a stretching exponent which 
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FIG. 1: T dependence of the diffusivity D, the viscosity n, 
and the product Drj/T for the cases iV m ax = 3, = 0.20 (top 
panel), and iV max = 4, <f> = 0.30 (bottom panel). Dotted lines 
correspond to the expected value (37rcr) _1 from the Stokes- 
Einstein relation. Dashed lines are fits to Arrhenius laws. An 
error bar is included for the viscosity at high T. 



takes values < f3 < 1. Experimental evidence for 
a collection of chemically and structurally very differ- 
ent glass-forming liquids Q shows that the smaller the 
fragility index (i.e., the closer the system is to strictly 
strong behavior), the closer to unity the values of (3 
are. Figure [3 shows that this is indeed the case for the 
iVmax model. The long-time dependence of the normal- 
ized coherent intermediate scattering function, 4> q {t) = 
{Pq{t)p~q(0))/{p q (0)p-q(0)), where p q {t) = E 2 =iexp[iq- 
ri(t)], can be well described by KWW fits with values of 
(3 > 0.85 in all the q range and for all studied <f>. Such 
a behavior is observed at all T where the system shows 
Arrhenius behavior. As shown in Ref. Q , these j3 values 
are very different from the ones typical of fragile liquids 
(/3 -0.5). 

Results reported in Figs. H an d HI provide convincing 
evidence that Eqs. JT][2Jl define a simple and satisfactory 
minimal model for a strong glass-forming liquid. 



FIG. 2: Symbols: Coherent intermediate scattering function 
for A^max = 4, (j> = 0.30, and T = 0.1, for different values of 
the wavevector q. Dashed lines are KWW fits. The stretching 
exponents (3 are indicated for the corresponding q's. 



b) Energy landscape 

Fig. [3]shows the T dependence of the potential energy 
(i.e., the energy of the typical bonding pattern) per parti- 
cle, E. In the present model, the potential energy of each 
configuration coincides with the energy of the bonding 
pattern and can be directly associated, in the Stillinger- 
Weber formalism, with the IS energy. Within the times 
accessible by the simulations, equilibrium states can be 
reached for configurations characterized by a number of 
broken bonds smaller than a 2 %, i.e. the lowest energy 
state, -Eft,, is approached from equilibrium simulations. 

From the low T behavior of E, one sees that the ap- 
proach to _Efb is well described by an Arrhenius law: 



E-Efb = ^ f 00 exp(-e f /fc B T). 



(15) 



The activation energy ef, determined by a free fit, is 
close to Mo/2- Indeed, by forcing the Arrhenius acti- 
vation energy to be exactly uq/2 a satisfactory repre- 
sentation of the data is recovered, with one simple fit- 
ting parameter Aoo . Figure [3] shows the result of a fit 
to E — Eft, = Aoo exp(— uo/2/cbT). The corresponding 
best-fitting values, with two (A^jUo) and one (^4oo) free 
parameters are reported in Table ||| The observed value 
uq/2 is consistent with theoretical predictions based on 
the thermodynamic perturbation theory developed by 
Wertheim [63| to study association in simple liquids. It 
is not a coincidence since in Wertheim's theory bonds are 
also geometrically uncorrelated. Similar values are also 
predicted by more intuitive recent approaches [64|. 

The clear low T Arrhenius dependence and the explicit 
value of the activation energy provide a convenient way 
to evaluate the T dependence of the energy for lower 
T. While in T, it might appear as a wide extrapolation 
procedure, we recall that in E the interpolation extends 
only over a small (2%) range of energies, between the fully 
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FIG. 3: T dependence of the potential energy per particle E 
for the cases iV max = 3, 4> = 0.20 (top panel), iV max = 4, 
<ji = 0.30 (middle panel), and -/V ma x = 5, <f> = 0.35 (bottom 
panel) . Full lines at low T are fits to Arrhenius behavior with 
activation energy uo/2 (see text and Table QJ. Dashed lines 
at high T correspond to linear behavior in T _1 (see text). 



bonded state (E — E^) and the lowest equilibrated state 
studied in simulations. Eq. (|15|l provides a convenient 
expression for the low T behavior of E and, by using Eq. 
lfT5)l. a way of calculating the total excess entropy down 
to the fully connected state. 

We note on passing that at intermediate temperatures 
the T dependence of E is consistent with a 1/T law, as 



(iVmax,0) 








(3, 0.20) 


2.69 


2.70 


0.506 


(3, 0.30) 


1.74 


1.72 


0.499 


(3, 0.35) 


1.38 


1.37 


0.498 


(4, 0.30) 


2.59 


2.73 


0.508 


(4, 0.35) 


2.15 


2.12 


0.498 


(4, 0.40) 


1.70 


1.67 


0.498 


(4, 0.45) 


1.28 


1.23 


0.495 


(4, 0.50) 


0.900 


0.884 


0.498 


(5, 0.35) 


3.14 


3.73 


0.539 



TABLE I: Fit parameters for the low T Arrhenius dependence 
of the potential energy E (see text). 
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FIG. 4: Full lines: T dependence of the isochoric configura- 
tional specific heat for several values of (N m&K , <f>). Dashed 
lines are an extrapolation to high T of the low T behavior 
A^eik^T' 2 exp(-e f /fc B T) (see text). 

expected for a Gaussian distribution of energy levels. The 
1/T law crosses to the Arrhenius dependence on cooling. 
This crossing has been also observed in the study of the 
T dependence of the IS energy in a realistic (atomistic) 
model for silica 0] . 

The low T Arrhenius dependence of E [Eq. 115|) ] has a 
practical implication in the T dependence of the isochoric 
configurational specific heat C£? n£ (T) = {dE/dT) v . 
Hence, from Eq. at low T we have C£? nf (T) = 

A^et fcg 1 T~ 2 exp(— et/k^T), which has a maximum at 
T = ef/2 w uo/4. Fig. 21 shows that indeed, numeri- 
cal data for Cy onf display a peak at T w 0.25. A peak 
in Cy nl (T) has also been observed in recent simulations 
of atomistic models of two different network-forming liq- 
uids: silica [HI and BeF 2 01 . 

We note that a strong correlation is observed between 
the T dependence of the diffusivity and the T dependence 
of the potential energy. Fig. shows, for N max = 4 
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FIG. 5: T dependence of the diffusivity (top panel) and 
the isochoric configurational specific heat (bottom panel) for 
A'max = 4, at several values of <f>. Dashed lines in top panel 
are Arrhenius fits. 



and several <j) values, that on cooling, D crosses to an 
Arrhenius law at T w 0.25. This temperature is the 
same at which the specific heat shows a maximum. This 
correlation holds at all the investigated cf> values. 

The quality of the low T Arrhenius fits for the diffu- 
sivities is worse at high tj>. Indeed, low T data at high 
4> show some bending (Fig. [SJl). This result suggests 
that the system will become more fragile with increas- 
ing density. This is not surprising, since the influence of 
the square well will be weaker at higher packing, and the 
system will approach a dense hard sphere liquid, which 
is a fragile system. 

Next we turn to the evaluation of the statistical prop- 
erties of the PEL, and more precisely the total excess 
entropy and its vibrational and configurational contribu- 
tions. As explained in Section III [Eq. (| 1 31) ] . the total 
excess entropy 5£* t is evaluated by isochoric integration 
from the hard sphere limit at a sufficiently high T IC { . We 
use a value T re f = 100. Fig. EJshows S°* t as a function of 
T and E, for different values of 4> and N max . Due to the 
presence of the interaction potential Q |2J , Sfot * s nega- 
tive (see also below). As expected, S^ t (T) decays to a 




-5.5 

, ex 
'tot 

-6 
-6.5 

-5 

-.ex -6 
tot 

-7 



o 


S»(E) 


• 


s ex (T) 




-2.5 



0.05 0.1 0.15 T 0.2 

-2.4 -2.3 -2.2 -2.1 -2 -1.9 E -1.1 



1 


1 i 1 


I 1 




O S«(E) 






• s»ro 








N =5 $ = 0.35 ~ 






max T 






i.i.i.i 



0.05 



0.1 



0.15 



0.2 0.25 T 0.3 



FIG. 6: E and T dependence of the total excess entropy over 
the ideal gas value, Stot, for the cases iV m ax = 3, (j) = 0.20 (top 
panel), N max = 4, cf> = 0.30 (middle panel), and N max = 5, 
cj> — 0.35 (bottom panel). Continuous lines for St£ t (E) are 
obtained as the sum of the fit functions for S^ b (E) [Eq. 1161 1 
and S con {(E) [Eq. 1191 1. Continuous lines for S^Jt(T) are 
parametrically obtained from the T dependence of E. 



constant value at low T, since the system is very close to 
the fully bonded state and no further structural changes 
are expected to occur. A glass transition temperature T g 
can be operationally defined as the T at which relaxation 
becomes longer than the simulation time. The fact that 
the system is so close to its lowest energy state already 
above this operational T g implies that a very small drop 
of the specific heat is expected at T g , consistently with 
experimental evidence in strong liquids |5j. 

The evaluation of the vibrational contribution S°f h [Eq. 
ifTTJl] requires the calculation of the integral over the cou- 
pling constants A. Fig. [7| shows the calculated A depen- 
dence of ^(J2iLi( r i — r° i ) 2 )\/N at several T for one spe- 
cific value of N max and <f>. Data for different iV max and 
4> display a similar behavior. At values A < A m j n « 10~ 6 
contributions to the integral in Eq. I|ll(l are negligible, 
and Amin is taken as lower-cut for integration. At large A 
values, A(^^ 1 (r i — r ) 2 )^/^ approaches the theoretical 
limit iksT/2. This value is reached at A max > 10 6 . Note 
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FIG. 7: A dependence of /3\(£*Li ( r * " ^) 2 )x/N, for iV max = 
4 and <f> = 0.30, at different temperatures. The horizontal line 
indicates the expected value 3/2 for harmonic behavior. 



that A = 10 6 corresponds to an average displacement per 
particle of the order of 10~ 4 for this T range. Hence 
the harmonic perturbation localize the particles in a well 
much narrower than the square well width A, so that 
the presence of the unperturbed potential is irrelevant 
for this upper A value. 

Fig. |H1 shows the T and E dependence of S°f h for the 
same values of <j) and N max of Fig. El Close to the fully 
connected state S^f h can be well described by a linear 
dependence on the energy (i.e. on the number of bonds): 
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S\ib{E) — SvTb(Efb) + lvib{E ~ Efb). 



(16) 



FIG. 8: As Fig. E]for the excess vibrational entropy, S^f b . 
Dashed lines for S^f h (E) are linear fits [Eq. I|16|l ]. Continu- 
ous lines for S^,f b (T) are parametrically obtained from the T 
dependence of E (see text). 



Interestingly, this linear dependence on the basin depth 
has also been observed in prev iously investigated models 
of supercooled liquids 0,13 El EE HH • 

Finally, the configurational entropy S con f is determined 
as Sfo t — S®f h . Fig. |5]shows S'conf for the same <j> and iV max 
of Figs. El and |H1 Some considerations are in order: as for 
the total entropy, S'conf approaches a constant value at 
low T. Interestingly enough, this constant value is signif- 
icantly different from zero. The fully connected network 
is thus characterized by an extensive number of distinct 
bond configurations ~ exp(-/VS CO nf )■ These different net- 
work configurations arise from different bond topologies, 
i.e., disorder is associated to the presence of closed loops 
of different number of bonds. 

To derive a functional form for the E dependence 
of S'conf, we start from the thermodynamic relation 
(dS/dE)v = 1/T, which in the present case can be writ- 
ten as: 



<9(S con f + S v ib) 



1 

T' 



d(E - Eo,) 

At low T, from Eq. JTSJ we obtain 1/T = -(l/e f ) \n[(E 



(17) 



Efbj/A^] and, making use of Eq. (|l(i|> we find: 

<9S C onf 1 , E — Epo 



d(E - Efb) 



- 7 vib--ln- 

ef a x 



(18) 



which, after integration, provides the E dependence of 
the configurational entropy: 



E — E\ 



fb 



Cf 



Sconf(T) — S con f(£'fb) 

ln[2(T - Efb)] + (E-Eb,), (19) 



where the constant 7 con f is given by: 

i H^aU 

7conf = 7vib H • (20) 

Cf €f 

As mentioned above, a satisfactory description of the low 
T Arrhenius dependence of E is provided by forcing a 
value uq/2 for the activation energy ef. Hence, we can 
make the changes ef — ► uo/2 and — > A^ in Eqs. (jT§l 
120(1 . These changes allow us to obtain a simple expression 
of Sconf in terms of the number of broken bonds. Hence, 
since E — Efb = NbhUo, we find: 

Sconf (E) = Sconf -2iV bh ln(2JV hb ) +7conf^bb- (21) 
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Seoul' Dashed lines for S cou {(E) are fits to Eq. 1191 . Contin- 
uous lines for S C onf(T) are parametrically obtained from the 
T dependence of E (see text). A typical error bar is shown in 
down panel. 



This expression suggests that the low-E,T dependence 
of the configurational entropy is controlled by a combi- 
natorial factor related to the number of broken bonds 
randomly distributed along the network. The derivation 
also shows how intimately the logarithmic dependence of 
the entropy is connected to the Arrhenius dependence of 
E at low T. The final expression for S con f(E) is very 
different from the quadratic energy dependence resulting 
from the Gaussian distribution of IS energies observed in 
models of fragile liquids HI |H H3 . 

To provide a further and analysis-free confirmation of 
the crossover at low T toward combinatorial statistics 
of the bonding energy states, we show in Fig. the 
T dependence of E for N max = 4 at several <j> values. 
We show a representation of E — Epo, both in linear and 
logarithmic scale, as a function of 1/T. We note that, for 
all </>, when E—Efc « 0.75 the 1/T law, which is expected 
to hold in a Gaussian landscape, breaks down. Similarly, 
when E — £fb « 0.3 the Arrhenius law sets in. The 
fact that both crossovers are, at most, weakly dependent 
on suggests that they are esentially controlled by the 
bonding statistics and that between E — Er, « 0.75 and 
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FIG. 10: T dependence of E — En, at different regions of the 
energy landscape. Top panel: Full lines correspond to behav- 
ior linear in 1/T. The horizontal dashed line indicates the 
departure of such behavior. Bottom panel: Full lines cor- 
respond to Arrhenius behavior. The horizontal dashed lines 
indicate the limits of Arrhenius and 1/T behavior. 



E — Eft, « 0.3 a crossover from Gaussian to logarithmic 
statistics occurs. 

Fig. El shows the E and T dependence of S CO nf for 
different values of N max . Eq. IjlOJI provides a good de- 
scription of the data. No fit parameters are involved in 
comparing the numerical estimates of S'conf and the pre- 
dictions of Eq. Ijl 90 > except for the constant S'conf (T"fb)- 
Note that 7 CO nf is not a fit parameter, but a function [Eq. 
iPty] of parameters defining E(T) and S°f h (E). Before 
discussing the calculated values for S con {(Efb), we present 
in Fig. ^2 the E dependence of the configurational and 
excess vibrational entropies for iV max = 4 and different 
packing fractions 4>. For all 4> 1 the same functional forms 
of Eqs. (|16|l and Ijl9|l are recovered respectively for S°f b 
and S'conf. The two quantities show an opposite trend. 
Curves for S con {(E) tend to collapse at high (f>, while 
those for S°f h (E) tend to collapse at low <j>. Interestingly, 
the configurational entropy just shifts with varying <p. 

Table [H] summarizes the results of the fits of the vi- 
brational and configurational entropies to Eqs. I|10|) and 
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FIG. 11: E dependence of S^f b and Sconf f° r JVmax = 4 and 
different values of <f>. Dashed lines in top and bottom panels 
are, respectively, fits to Eqs. 1161 and I19H . 



Ijl9|l for the studied range of control parameters. Data 
shown in the table help discussing the N max and <f> de- 
pendence of the entropy of the fully bonded state. A 
trend in the direction of increasing S C onf(-Efb) on increas- 
ing iV max is observed very clearly in the comparison be- 



tween N n 



3 and N„ 



4. Much weaker is the 



trend between N max — 4 and iV max = 5. A similar 
weak trend is observed for N max = 4 on increasing <f>. 
The weak increase of SconK-Efb) with increasing <fi sug- 
gests that when the system is compressed, neighboring 
particles progressively enter in the interaction range of 
a given one, yielding a major variety of local configu- 
rations of the bonding pattern, and consequently, more 
topologically distinct fully bonded networks. The trend 
of the excess vibrational entropy suggests that the avail- 
able free volume for a given bonding pattern decreases 
with increasing <f>. 

A summary of the landscape analysis for all studied 
<fi is shown in Fig. ^]for N max = 4. The figure shows 
the 4> dependence of Sf* t , S*°f b and S^onf for several dif- 
ferent low T isotherms, all in the Arrhenius region of 
energies. It clearly emerges that the significant reduc- 
tion of S^ t on increasing <j> arises essentially from the vi- 
brational component. We also note that the HS relative 



TABLE II: Parameters denning the configurational [Ea. 119^ ] 
and excess vibrational [Eg. 1161 1 entropy for the studied values 
of iVmax and <f>. 



contribution to S°* t increases on increasing 



according to Eq. 



Hence, 

(JTJ, Sfxty = 0.30) = -1.90, while 
0.50) = —5.0. By comparing with data in Fig. 
1121 it is clear that S%* t is dominated by the square-well 
and hard-sphere contributions at respectively low and 
high packing fraction. 

It is also interesting to observe that, within the preci- 
sion of the data, S con { shows a weak maximum, shifting 
to higher <fi on cooling. To test wether the presence of a 
maximum of .S^onf has some effect on the dynamics we 
also show in Fig. ^|the behavior of the diffusivity along 
the same isotherms. We note that D is monotonic in <fi 
and hence that the maximum in the configurational en- 
tropy does not provoke a maximum in the diffusivity. We 
also note that an isochoric plot (not shown) of log D vs. 
[T^conf (T)] _1 provides a satisfactory linearization of the 
data, as suggested by the Adam-Gibbs theory 66] . This 
is not inconsistent with the observed Arrhenius depen- 
dence of D, since the T dependence of S con f(T) is only 

at mOSt 20% Of Sconf (-Eft,). 

V. CONCLUSIONS 

This article reports an explicit numerical calculation 
of the potential energy landscape for a simple model of 
strong liquids. It shows that it is possible to calculate 
with arbitrary precision the statistical properties of the 
landscape relevant to the behavior of the system at low 
T, when all particles are connected by bonds. The model 
can be seen as a zero-th order model for network-forming 
liquids, capturing the limited valency of the interaction 
and the open structure of the liquid. By construction, it 
misses all geometric correlations between different bonds 
which are present in network-forming materials. The 
simplicity of the model has several advantages, some of 
which of fundamental importance for an exact evalua- 
tion of the landscape properties. Hence, since angular 
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constraints between bonds are missing, it is possible to 
equilibrate the system to very low T, reaching configu- 
rations which are essentially fully bonded. At the lowest 
studied T, less than 2% of the bonds are broken in av- 
erage. Differently from other studied models with fixed 
bonding site geometries j^, El> EH > the absence of geo- 
metric constraints makes it possible to reach almost fully 
bonded states in a wide range of densities. Moreover, the 
use of a square well interaction as bonding potential has 
the advantage that the energy of the fully bonded state 
(the lowest possible energy) is known. 

The present model neglects completely interactions be- 
tween particles which are not nearest neighbors. The en- 
ergy of a particle is indeed fully controlled by the bonds 
with the nearby particles. This element of the model fa- 
vors a sharp definition of the energy and a clear-cut def- 
inition of basins. On the other hand, in network-forming 
liquids, bonding is often of electrostatic origin and inter- 
actions are not limited to the first shell of neighbors. This 
produces a much wider variety of local environments and, 
as a consequence, a spreading of the energy levels. These 
residual interactions, if smaller than the bonding inter- 
actions, will only contribute to spread the distribution 
of energy states without changing the landscape features 
(down to T of the order of the energy spreading) . 

The use of square- well interactions has a major advan- 
tage in relation to the possibility of precisely calculating 
landscape properties, since the energy of the system be- 
comes a measure of the number of bonds. A basin in con- 
figuration space can be identified as a bonding pattern, 
and a transition between different basins becomes asso- 
ciated to bond forming and breaking. Under these con- 
ditions, the basin boundaries are properly defined. We 
have shown that the method of the Perturbed Hamilto- 
nian can be extended to the present model, providing a 
formally exact method to evaluate the vibrational com- 
ponent of the free energy. This is a relevant achievement, 
since the evaluation of the vibrational entropy is the only 
weak point in all estimates of landscape properties in 
models with continuous potentials Indeed, when 

the potential is continuous, the constraint of exploring a 
fixed basin can not be implemented unambiguously in 
the Frenkel-Ladd method due to the difficulty of detect- 
ing crossing between different basins. Such a difficulty is 
not present in the square-well potential, since crossing of 
a basin is detected by a finite energy change. 

The two relevant features observed in this study are: 
(i) A residual value of the configurational entropy for 
T —> 0, associated to the exponentially large number of 
distinct fully connected bonding patterns, (ii) A loga- 
rithmic dependence of the number of bonding patterns 
on the number of broken bonds [Eq. J2J]. These two 
features are common to all investigated values of N max 
and to all studied 4>. A consequence of the logarithmic 
landscape statistics is the absence of a finite temperature 
at with the lowest energy state is reached [7(| . Indeed 
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from Eq. (jT§|l it is found that dS/dE\E=E fb = oo, and 
hence is reached only at T = 0. 

According to the picture emerging from this study, 
strong liquid behavior is connected to the existence of an 
energy scale, provided by the bond energy, which is dis- 
crete and dominant as compared to the energetic contri- 
butions coming from non-bonded next-nearest neighbor 
interactions |71j . It is also intimately connected to the ex- 
istence of a significantly degenerate lowest energy state, 
favoring the formation of highly bonded states which can 
still entropically rearrange to form different bonding pat- 
terns with the same energy. Of course, the specific value 
of /Sconf (Eft,) in the A max model provides an upper bound 
to the value expected in network- forming liquids [zSIz^l, 
since the absence of angular constraints significantly in- 
creases the number of geometric arrangements of the par- 
ticles compatible with a fully bonded state. Recent es- 
timates in glassy water ,741, which forms a tetrahedral 
disordered network, suggest a residual value of the con- 
figurational entropy of the order of ~ 0.3/cb ■ Hence the 
localization of the bonding sites at specific locations, and 
the associated geometrical correlations, do produce a sig- 
nificant reduction of .Sccmf (£fb)- 

Results reported in this article suggest that strong and 
fragile liquids are characterized by significant differences 
in their potential energy landscape properties. A non- 
degenerate disordered lowest energy state and Gaussian 
statistics characterize fragile liquids, while a degenerate 
disordered lowest energy state and logarithmic statistics 
are associated with strong liquids. These results ratio- 
nalize previous landscape analysis of realistic models of 
network- forming liquids |3lj . and the recent observation 
by Heuer and coworkers that the breakdown of Gaussian 
landscape statistics is associated with the formation of a 
connected network [43j. While in atomistic models the 
lowest energy state is not known, and the very long equi- 
libration times prevent an unambiguous determination of 
its degeneracy, both these quantities are accessible in the 
present simple model. 

A last remark concerns the limit of the A max value 
for which the fully connected state can be reached. 
The possibility of approaching the fully bonded state is 
limited by the possibility of avoiding the (T — (f) region 
where liquid-gas phase separation is present. It has 
recently been observed [13)113, E3 that the region of un- 
stable states expands on increasing A max and essentially 
covers, at low T, the entire accessible <p range when 
A max > 6. For these large A max values, slowing down 
of the dynamics is observed only at very large <j) and it 
is essentially controlled by packing considerations, not 
by bonding. In this respect, bond-controlled dynamics 
are observable only when the valence of the interparticle 
interaction is limited. 
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